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Abstract. Transient simulations of a resonant tunneling diode oscillator are presented. 
The semiconductor model for the diode consists of a set of time-dependent Schrodinger 
equations coupled to the Poisson equation for the electric potential. The one-dimensional 
Schrodinger equations are discretized by the finite-difference Crank-Nicolson scheme using 
memory-type transparent boundary conditions which model the injection of electrons from 
the reservoirs. This scheme is unconditionally stable and reflection-free at the boundary. 
An efficient recursive algorithm due to Arnold, Ehrhardt, and Sofronov is used to im- 
plement the transparent boundary conditions, enabling simulations which involve a very 
large number of time steps. Special care has been taken to provide a discretization of 
the boundary data which is completely compatible with the underlying finite-difference 
scheme. The transient regime between two stationary states and the self-oscillatory be- 
havior of an oscillator circuit, containing a resonant tunneling diode, is simulated for the 
first time. 



I. Introduction 

The resonant tunneling diode has a wide variety of applications as a high-frequency and 
low-consumption oscillator or switch. The resonant tunneling structure is usually treated as 
an open quantum system with two large reservoirs and an active region containing a double- 
barrier heterostructure. Accurate time-dependent simulations are of great importance 
to develop efficient and reliable quantum devices and to reduce their development time 
and cost. There exist several approaches in the literature to model a resonant tunneling 
diode. The simplest approach is to replace the diode by an equivalent circuit containing 
nonlinear current- voltage characteristics [28] . Another approximation is to employ the 
Wannier envelope function development Other physics-based approaches rely on the 
Wigner equation [T2J EZ], the nonequilibrium Green's function theory [TJ)J [22], quantum 
hydrodynamic models [191 ETJ EI], and the Schrodinger equation (TT], [131 EH 132] . 

In this paper, we adopt the latter approach and simulate the time-dependent behavior of 
a resonant tunneling diode using the Schrodinger-Poisson system in one space dimension. 
In this setting, the electrons are assumed to be in a mixed state with Fermi-Dirac statistics 
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and the electrostatic interaction is taken into account at the Hartree level. Each state 
is determined as the solution of the transient Schrodinger equation with nonhomogeneous 
transparent boundary conditions. The Schrodinger equations are discretized by the Crank- 
Nicolson finite-difference scheme and coupled self-consistently to the Poisson equation. The 
main originality of this paper is the adaption (and slight improvement) of existing numerical 
techniques from, e.g., [U El [31], to a long-time simulation of a high-frequency oscillator 
circuit containing a resonant tunneling diode. Our changes in the techniques are necessary 
to achieve simulations without spurious oscillations in the numerical transient solution. In 
the following, we detail the techniques used as well as the novel features. 

First, we consider the one-dimensional stationary problem, since it builds the basis for 
the transient simulations. The stationary transparent boundary conditions are discretized 
in such a way that their discrete version is compatible with the underlying finite difference 
discretization, as proposed in [3]. Thereby, any (numerical) spurious oscillation is elim- 
inated, which would otherwise propagate in the transient simulations. In the literature 
|10| [31] . a modified version of the potential energy is used to overcome problems of nu- 
merical convergence. Physically interpreted, this model introduces artificial surface charge 
densities at the junction interfaces of the tunneling diode. We are able to solve the original 
problem. This represents an improvement compared to the simulations in [3T], where the 
modified model is also employed for the time-dependent case. 

Second, the time-dependent Schrodinger-Poisson system with transparent boundary con- 
ditions is approximated. Since these boundary conditions are of memory type [U [8] , their 
numerical implementation requires to store (and to use) the boundary data for all the 
past history. For this reason, simulations involving longer time scales are extremely costly. 
This explains why simulations in the literature [TOl IT3] [31] have been restricted to some 
picoseconds only. We solve this problem by using a fast evaluation of the discrete convo- 
lution kernel of sum-of-exponentials, which has been presented in [B] and employed in [3] 
on circular domains. To our knowledge, this rather new numerical technique has not been 
applied to realistic device simulations so far. 

A challenge in the transient simulations results from the large number of wave func- 
tions which need to be propagated, accounting for the energy distribution of the incoming 
electrons. Each state is provided with transparent boundary conditions, which raises the 
computational cost sharply. To cope with a large number of Schrodinger equations to 
be solved, we developed a parallel version of our solver utilizing multiple cores on shared 
memory processors. This enables us to present, for the first time, simulations to the 
Schrodinger-Poisson system for large times up to 100 ps (ps = picosecond) with reasonable 
computational effort (compared to 5ps in [21], 6ps in [T3], and 8ps in [TP]). 

Another novelty in this paper concerns the discretization of nonhomogeneous discrete 
transparent boundary conditions. They are necessary to describe continuously varying 
applied potentials (as in an oscillator circuit). It is well known that, using a suitable 
gauge change, one can get rid of the transient potential [2]. Corresponding nonhomoge- 
neous transparent boundary conditions can be found in [10J. In numerical simulations, 
however, we observed that these boundary conditions may lead to unphysical distortions 
in the conduction current density. The reason is that the considered discretization of the 
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gauge change is not compatible with the underlying finite-difference scheme. Therefore, 
we suggest a new discretization which is derived from the Crank-Nicolson time integra- 
tion scheme. Our approach completely removes these numerical artifacts, and we show 
that the total current density is now perfectly conserved. We stress the fact that our dis- 
cretization is completely consistent with the underlying Crank-Nicolson scheme inheriting 
its conservation and stability properties. 

Third, the numerical results allow us to identify plasma oscillations in a certain time 
regime of the resonant tunneling diode and to estimate the life time of the resonant state. 
We present, for the first time, simulations of a high-frequency oscillator circuit containing 
a reconant tunnelding diode, based on a full Schrodinger-Poisson solver with transparent 
boundary conditions. Simplified tunneling diode oscillators have been considered in [2H1 
1291 30J. Our approach enables us to observe the complex spatio-temporal behavior of 
macroscopic quantities inside the resonant tunneling diode in an unprecedented way. 

The paper is organized as follows. In Section [2j we detail the algorithm of the station- 
ary problem. The transient algorithm is described in Section |3j In Section [4] we consider 
numerical experiments for constant applied voltage and time-dependent applied voltage. 
Furthermore, the numerical convergence related to the approximation of the discrete con- 
volution kernel by sum-of-exponentials is investigated. Finally, we present high-frequency 
oscillator circuit simulations in Section |U 

2. Stationary simulations 

The steady state is the basis for the transient simulations. Therefore, we discuss first 
the stationary regime. 

2.1. Schrodinger-Poisson model. We assume that the one-dimensional device in (0, L) 
is connected to the semi-infinite leads (— oo, 0] and [L, oo). The leads are assumed to be in 
thermal equilibrium and at constant potential. At the contacts, electrons are injected with 
some given profile. We suppose that the charge transport is ballistic and that the electron 
wave functions evolve independently from each other. The one- dimensional device consists 
of three regions: two highly doped regions, [0, aj and [a 6 , L], with the doping concentration 
n D and a lowly doped region, [ax, a 6 ], with the doping density n 2 D (see Figure]!]). The middle 
interval contains a double barrier, described by the barrier potential 




V 




for x G [a 2 , a 3 ] U [a 4 , a 5 ], 
else. 



The doping profile no is defined by 




else. 



for x G 




The parameters are taken from [TP] I3"T] : 



ax = 50 nm, a 2 = 60 nm, 
a 4 = 70 nm, a 5 = 75 nm, 
L = 135 nm, n l D = 10 24 m 
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03 = 65 nm, 
a$ = 85 nm, 
n 2 n = 5 • 10 21 
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and the barrier height is Vq = 0.3eV. 
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Figure 1. Barrier potential and doping profile of a double-barrier het- 
erostructure. 



The Coulomb interaction is taken into account at the Hartree level, i.e. by an infinite 
number of Schrodinger equations 



(1) 



h 2 d\ 



(x) + V(x)(j> k {x) = E{k)<f> h {x) 



x G 



2m dx 2 

self-consistently coupled to the Poisson equation, 

(2) dx 2 " e 

Keif(O) = 0, V seli (L) 

where V = Vbarr + Keif is the potential energy. The physical parameters are the reduced 
Planck constant h, the effective electron mass m, the elementary charge e, and the permit- 
tivity e = e r eo, being the product of the relative permittivity e r and the electric constant 
Eq. Furthermore, U > denotes the applied voltage, and the electron density is defined by 



n D ), x e (0,L), 
-eU, 



(3) 



n [Keif] {x) 



g(k)\Mx)\ 2 dk. 



The injection profile g(k) is given according to Fermi-Dirac statistics by 



(4) 



g(k) 



mk B T 
2n 2 h 2 



In 1 + exp 



Ep - h 2 k 2 /(2m) 
k^fo 



where ks is the Boltzmann constant, To the temperature of the semiconductor and Ep the 
Fermi energy (relative to the conduction band edge). In all subsequent simulations, we 
use, as in [31], e r = 11.44, T = 300 K, E F = 6.7097 • 10 -21 J, and the effective mass of 
Gallium arsenide, m = 0.067m e , with m e being the electron mass at rest. 
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In order to define the total electron energy E(k) depending on the wave number k G M, 
we need to distinguish the cases k > and k < 0. For k > 0, the electrons enter from the 
left, and we have E{k) = H 2 k 2 /(2m). The wave function in the leads is given by 

fc (x) = e lkx + r(k)e~ ikx , x < 0, 



Ma; J 



t(A;) exp (iy/2m{E{k) -V(L))/h 2 xj , x > L. 



Eliminating the transmission and reflection coefficients t(k) and r(k), respectively, the 
boundary conditions 



(5) 0' fc (O) + ikMQ) = 2ik > <t>'k{L) = i^2m{E{k) - V(L))/h 2 ML) 

follow. For k < 0, the electrons enter from the right. The total energy is given by 
E(k) = h 2 k 2 /(2m) — ell, and the wave function in the leads reads as 

<f>k( x ) — t(k) exp y—iy/ 2mE(k) / h 2 xj , x < 0, 
<j) k (x) =e ikx + r(k)e~ ikx , x > L. 
This yields the boundary conditions 

(6) 0' fc (O) = -i^2mE(k)/H 2 (j) k (0), (f>' k (L) + ik<f> k (L) = 2ike ikL . 

Summarizing, the stationary problem consists in the Schrodinger equation ([T]) with 
the transparent boundary conditions ^ coupled to the Poisson equation ^ via the 
electron density We remark that the existence and uniqueness of solutions to a 

Schrodinger-Poisson boundary- value problem similar to ([l])-([6]) has been shown in [9]. 

2.2. Discrete transparent boundary conditions. We recall the finite-difference dis- 
cretization of the stationary Schrodinger equation with transparent boundary conditions 
[I]. Using standard second-order finite differences on the equidistant grid Xj = jAx, 
j G {0, J}, with xj = L and Ax > 0, we find for the grid points located in the 
computational domain, 

2m(Ax) 2 

(7) <f> j+1 - 2^ + + K h2 } (E(k) - Vj)^ = 0. 

It is well known that a standard centered finite-difference discretization of the bound- 
ary conditions ^ and ^ may lead to spurious oscillations in the numerical solution 
In principle, the numerical errors can be made as small as desired by choosing Ax suffi- 
ciently small. However, since the stationary solutions will serve as intitial states in our 
transient simulations, we need to avoid any spurious oscillation, which would otherwise be 
propagated with every time step. 

For this, we apply (stationary) discrete transparent boundary conditions compatible with 
the finite-difference discretization ([7]) as proposed in [3]. For the sake of completeness, we 
review the derivation. Note that the final discretization is equivalent to the discretization 
([7]) extended to the whole space, i.e. for j G Z. 
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In the semi-infinite leads j < and j > J, the potential energy is assumed to be 
constant, 

V = for j < 0, 

Vj = -eU for j > J. 

Then ^ reduces to a difference equation with constant coefficients which admits two 
solutions of the form <f>j = (a^jY, where 

± m(E(k)-V 0J )(Ax) 2 
«o,J = 1 jr 2 



+ j2m(E(k) - V ,j)(Ax) 2 m 2 (E(k) - H,j) 2 (Ax) 4 



h 2 h 4 

Here, E(k) — Vo,j corresponds to the kinetic energy Z?Qj(fc) in the left or right lead. In 
case #oj(fc) > 0, the solution is a discrete plane wave and (Ax) 2 < 2h 2 /(m(E(k) — Vo,j)) 
is needed to ensure |«o,j| = 1> which in practise is not a restriction. In case E^j(k) = 0, 
the solution is constant. Depending on the applied voltage, E^{k) might also become 
negative. In that case, the solution is decaying or growing exponentially fast and we select 
the decaying solution as it is the only physically reasonable solution. 

In practice, we start with the calculation of the total energy E(k) = E^™(k) + V^j. For 
electrons coming from the left contact we have E(k) = E^ in (k). As the incoming electron 
is represented by a discrete plane wave, E^ n (k) is positive but, depending on the applied 
voltage, Ej m (k) might be positive, zero or negative. For electrons coming from the right 
contact, we have E(k) = E ] f n {k) — eU. Again, the incoming wave function is a discrete 
plane wave, i.e., Ej m (k) > but nothing is said about E^ n (k). At this point it should be 
noted that the kinetic energy of the incoming electron needs to be computed according to 
the discrete dispersion relation 

E kkl (k) = —^—(l-cos(kAx)), 
m{L\x) z 

which follows after solving the centered finite-difference discretization of the free Schrodin- 
ger equation 

h 2 d 2 

-—^-e lkx = E kin (k)e ikx . 
2m ax 2 

In the limit Ax — > 0, we recover the continuous dispersion relation E km (k) = h 2 k 2 /(2m). 

Let us consider a wave function entering the device from the left contact (k > 0). For 
j < 0, the solution to (j7l) is a superposition of an incoming and a reflected discrete plane 
wave, <pj = ft + B(3~i , where (3 = a . We eliminate B from = + B(3, 4> = 1 + B 
to find the discrete transparent boundary condition at Xo: 

-/rV-i + 0o = i-/r 2 - 

For j > L, the solution to ^ is given by <pj = CV with 7 = cej. This means that 
0j + i = C r y J+1 = 70 j, and the boundary condition at xj becomes 

0j - 7"Vj+i = 0. 



TRANSIENT SCHRODINGER-POISSON SIMULATIONS 



7 



Summarizing, we obtain the linear system Acj) = b with the tridiagonal matrix A consist- 
ing of the main diagonal (-/3 _1 , -2 + 2m(Ax) 2 (E(k) - V ) /ti 2 , ...,-2 + 2m(Ax) 2 (E(k) - 
Vj)/h 2 , — 7 _1 ) and the first off diagonals (1, . . . , 1). The vector of the unknowns is given 
by — • • • ^ J+i) T an d b represents the right-hand side b = (1 — (3~ 2 , 0, . . . , 0) T . 

The case of a wave function entering from the right contact (k < 0) works analogously. 

2.3. Solution of the Schrodinger-Poisson system. We explain our strategy to solve 
the coupled Schrodinger-Poisson system. To this end, we introduce the equidistant energy 
grid 

(8) K = {-k M , -k M + Ak, . . . , -Ak, +Ak, . . . , k M - Ak, k M }, K:=\K\. 
The electron density rt3| is approximated by 

n d isc[V sc u](x) = Ak^g(k)\(j) k {x)\ 2 , 

keK. 

where the Fermi-Dirac statistics g(k) is defined in Q and the functions 0^ are the scat- 
tering states, i.e., the solutions to the discretized stationary Schrodinger equation Q with 
discrete transparent boundary conditions as described in Section 22 This approximation 
is reasonable if Ak is sufficiently small and kj^ is sufficiently large. In the numerical simu- 
lations below, we choose K = 3000 and, as in PH Section 5], k M = ^2m(E F + 7k B T )/h, 
recalling that E F = 6.7097 • 10~ 21 J and T = 300 K. 

The discrete Schrodinger-Poisson system is iteratively solved as follows. We set V = 



r(p) where V {p) 
self,C/' WIlele v se li,i 

we compute a set of quasi eigenstates 



Vbarr + Keif ui wnere V^l v is the p-th iteration of V^ e if for the applied voltage U . Given V, 



}ke)c- This defines the discrete electron density 



\V^ u ] = Ak"£9{mi P \x)\ 2 . 



kdK. 



The Poisson equation is solved by employing a Gummel-type method 
d 2 



T/(P+1) 

dx 2 seU ' u 



e 
e 



nlV^n] exp 



v se\i,U v aeU,U 



yrcf 
V sclf 



n D 



(p+i) 



(0)=0, V™{L) 



self,f/ 



-eU. 



The idea of the Gummel method is to decouple the Schrodinger and Poisson equations 
but to formulate the Poisson equation in a nonlinear way, using the relation between 
the electron density and electric potential in thermal equilibrium. The parameter 
can be tuned to reduce the number of iterations; we found empirically that the choice 
Keif = 0.04 eV minimizes the iteration number. If the relative error in the £ 2 -norm is 
smaller than a fixed tolerance, 



(9) 



(p+i) 



self,!/ 



V. 



(p) 



scli,U 



(P+l) 
self,!/ 



< 5, 
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we accept V^if := and {0^, p+1 ^}fc e /c as the approximate self-consistent solution. Oth- 

erwise, we proceed with the iteration p + 1 — > p + 2 and compute a new set of scattering 
states. The procedure is repeated until ^ is fulfilled. We have choosen the tolerance 
5 = 1(T 6 . 

For zero applied voltage we use 0mV = OmV to start the iteration. Only 7 iterations 

are needed until criterion ^ is fulfilled. As a result we obtain V^j? 0mV , which is depicted 
in Figure [2] (left part, solid line). 

Numerical problems arise when non-equilibrium solutions are computed. As an example 
we consider the case of a small applied voltage U = 1 mV. To start the iteration process we 
use the previously computed solution, i.e., we set V^ {lmV := 0mV . The next iterations 



are illustrated in Figure [2] (left part, dashed lines). Obviously they do not converge and 
are physically not realistic. This phenomenon is a well-known in the literature \XT\ [T5] and 
is believed to be related to the absence of inelastic processes in the Schrodinger-Poisson 
equations. 

In the literature [101 |31], a modified version of the Schrodinger-Poisson equations is 
employed to overcome this problem. The modification concerns the description of the 
potential energy in the Poisson equation. For this, we write the Poisson equation ^ as 
follows: 

(M " m(0,L), V o (0)=0, V (L) = -ell, 



dx 2 
d 2 V x _ 
dx 2 e 



n-n D ) in(0,L), ^(0) = 0, V\(L) = 0, 



i.e., the self-consistent potential is V^ e if = Vo + Vi. The first boundary- value problem can be 
solved explicitly: Vq(x) = —eUx/L, x G [0, L\. In pill ED], the linearly decreasing potential 
Vq has been replaced by the ramp-like potential 

— fX — (2i \ 

(10) V (x) = -B l [aii06) + l [a6>oo) , x G [0, L], 

where 1/ is the characteristic function on the interval / C IRfsee Figure [l] for the definition 
of a\ and as). The function Vo + H.arr is illustrated in Figure |2| (right part, dotted line). The 
potential energy is then given by V — Vq + V\ + H.arr- Using this modified physical model, 
the above Gummel iteration scheme for the Poisson equation for V\ converges without any 
problems, see Figure [2] (right part, solid line), even for large applied voltages. However, we 
will see below that the results from the modified model differ considerably from the results 
obtained by the original Schrodinger-Poisson model. Furthermore, the potential energy is 
no longer differentiable at ai and clq. This may be interpreted as a model of surface charge 
densities at the interfaces which, however, are not intended in the model. 

In fact, we are able to solve numerically the original Schrodinger-Poisson problem. To 
this end, the applied voltage needs to be increased in small steps. We found that the 
starting potential in each step needs to be initialized carefully. More precisely, given the 
self-consistent solution V^ e if,u f° r the applied voltage U, we wish to compute a self-consistent 
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Figure 2. Left part: Solid line: self-consistent solution V^ e if for U = OmV, 
found after 7 iterations. Dashed lines: divergent approximations for U = 
lmV. Dotted line: barrier potential. Right part: Solid line: self-consistent 
solution V\ for U = 250 mV according to approximation (10). Dotted line: 
sum of the barrier potential and a ramp-like potential. 



solution with the applied voltage U + AU. In each step we choose 
(11) V$t,u + *n(z) ■= " At^^l [w] 

to start the iteration. For U = OmV and AU = 25 mV, the Gummel scheme converges 
to a physically reasonable solution after 7 iterations (i.e., (|9j is fulfilled). Some iterations 
are shown in Figure [3j We observed that a voltage step AU < 30 mV leads to convergent 
solutions also for large applied voltages. 

In order to compare the original Schrodinger-Poisson model with the model using ap- 



proximation (10), we computed the current- volt age characteristics shown in Figure 4 Here 



the (conduction) current density 
(12) J cond = ^-£g{k)Im \4>t^J dk 

is approximated by a simple quadrature formula using symmetric finite differences to com- 



pute d<fik/dx. Figure [4] shows that the results differ considerably, i.e., the choice (10) leads 
to different results than those computed from the original model. Therefore, we employ 
the original potential energy in the transient simulations in the next subsection. 

3. Transient simulations 

In this section, we detail the numerical discretization of the transient Schrodinger equa- 
tions 

(13) jft_|^ = __-^ + y(.,t)^ fcj ^(.,0) = ^, xe [o,L], t>0, ke)C, 
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Figure 3. Some iterations computed according to (11). 




0.3 0.4 
applied voltage in V 

Figure 4. Current- voltage characteristics. The solid line corresponds to our 
solution of the original stationary Schrodinger-Poisson system. The dashed 



0.6 



line is obtained with the modified model using approximation (10). 



with discrete transparent boundary conditions, where /C is defined in (J8J). To simplify the 
presentation, we skip in this section the index k. 
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3.1. Nonhomogeneous discrete transparent boundary conditions. The transient 



Schrodinger equation (13) is discretized by the commonly used Crank-Nicolson scheme: 



, (n+l) , (-T3 o i TA( n + 1 / 2 )^ /("+!) i /(™+l) 

y)_ x + (iR-2 + wv> )ip) + ip j+1 
(14) V j j j j 

= -iifX +(iR+2- wV} n+1/2) )^ n) - ^ 



(n) 
+1' 



where tpj approximates ip(xj,t n ) with Xj = jAx and t n = nAt (j G Z, n G N ), y^"^ 1 / 2 ) 
approximates V(jAo;, (n + 1/2) At), and i? = Am(Ax) 2 /(hAt), w = —2m(Ax) 2 /h 2 . Under 
the assumptions that the initial wave function is compactly supported in (0, L) and that 
the applied voltage vanishes, V(x, t) = for x < and x > L, t > 0, it is well known (see, 



e.g., [H E]) that transparent boundary conditions for the Schrodinger equation (13) read 
as 



d5a) w M= M e -^*rmA dT , 

y ' dx y ' V nh dt J Q 

d5b) ^ i , 1) ._.fe-^rp, T , 

y J dx y ' 1 V nh dt J 

The (homogeneous) discrete transparent boundary conditions, based on the above Crank- 
Nicolson scheme, are given as follows (see [3] for the derivation): 

n 

(16a) ^+ 1 )_ s (o)^+i) = ^ s (n+i-^)_^) ) n > , 

n 

(16b) ^)- s W^) = Y j S^ +l - l ^f-^ n>0, 

with the convolution coefficients 



(17) ^ = ^_ j _j4 M1+ ^i + j _j 4M + ae 2n _ i 

and the abbreviations 



arctan 4 V = , = « = - ^ 2 (i2 2 + 16) e^/ 2 . 



Here, P n denotes the nth-degree Legendre polynomial (P-% = P-i = 0), and 5 n j is the 



Kronecker symbol. In practice, the coefficients defined in (17) are computed with an effi- 
cient three-term recursion, relying on the three-term recursion of the Legendre polynomials 
[T5] . The Crank-Nicolson scheme along with these discrete boundary conditions yields an 
unconditionally stable discretization which is perfectly reflection-free [3j Hj. 

Next, let the initial wave function be a solution to the stationary Schrodinger equation 
with energy E and let the exterior potential at the right contact be given by a time- 
dependent function, V(x,t) = —eU(t) for x > L, t > 0. This leads to nonhomogeneous 
transparent boundary conditions [1] . We describe our strategy to discretize these boundary 
conditions. Our approach is motivated by that presented in [TUJ Appendix B], but we 
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suggest, similarly as in [1], a discretization of the gauge change which is compatible with 
the underlying finite-difference scheme. Additionally, our approach requires only a single 
set of convolution coefficients instead of two. 

First, we derive a nonhomogeneous discrete transparent boundary condition at Xj = 
L. To this end, we consider the difference between the unknown wave function ip and 
the time evolution of the scattering state, <f>{x) exp(iEt/h), in the right lead [L, oo). We 
employ a gauge change to get rid of the time-dependent potential Vi(t) = —eU(t). As a 
consequence, the function if)(x,t) exp(i CVi r (s)ds/K) solves the free transient Schrodinger 
equation in [L, oo). Using a similar gauge change, a straightforward computation shows 
that cj)(x) exp(— i(E — VL(0))t/fr) solves the free Schrodinger equation in [L, oo) as well. 
Hence, 



i /* \ ii 



(18) <p(x,t) :=V>(z,*)exp( - / V L (s)ds J - 0(x) exp ( --(E - V L (0))t 



o 

for x G [X, oo) solves the free transient Schrodinger equation. Furthermore, tp(x,0) = 
for all x G [L, oo). Therefore, we could apply (15b) to derive a nonhomogeneous trans- 
parent boundary condition at the right contact. Instead we replace tp{x,t) by some ap- 
proximation (pj and subsequently apply (16b) to derive a discrete boundary condition 
compatible with the Crank-Nicolson scheme. The question is how to approximate the 

& 



quantities exp(i $*Vi,(s)ds/K) and exp(— i{E — Vj,(0))t/^). Indeed, the ad-hoc discretiza- 
tion for t = n/St, 



exp (i fv L {s)ds\ ^exp (^V^At 
(19) v ' ' 'it / 



•' \ / i . ■ - r (o)> 



exp l--(E- V L (0))tJ = exp l--(E- V^^nAt 

where = Vi(£At), is not derived from the underlying finite-difference discretization, 
causing unphysical numerical reflections at the boundary. In principle, these reflections 
can be made arbitrarily small for At — > 0. However, for practical time step sizes, the 
calculation of the current density would be still distorted. Our (new) idea is to apply 
a Crank-Nicolson discretization to a differential equation satisfied by exp(i J^VL(s)ds/h). 
Indeed, this expression solves 

dp i 

-(t) = j_V L (t)e(t), e(0) = 1. 
The Crank-Nicolson discretization of this ordinary differential equation reads as 

£ (n+l) = £ (n) + A 4yW)( e (^l) + e (n)) s £ (0) = L 

This recursion relation can be solved explicitly yielding 

eW = exp (^i arctan ^^ +1/2) ^) V n G N . 
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2* arctan f— lf +1/2) 



l V ? +x ' 2) At + 0{{Atf) 



reveals that in the limit At — > 0, the ad-hoc discretization in (19) coincides with the 
discrete gauge change which is derived from the Crank-Nicolson time-integration method. 
Analogously, exp(— %{E — Vi,(0))t/h) needs to be replaced by 



(n) 

7j 



exp 



exp 



n-l 

2i arctan 



2in ^arctan 



A* 



E^jj exp \ arctan {^^L^j 



At, 

~2h 



(0) 



arctan 



At 

~2h 



-E 



n G N n - 



Thus, the discrete analogon of cp in definition (Jl8|) is given by 



Replacing by tpj l> in ( 16b), we obtain the desired nonhomogeneous discrete transparent 
boundary condition at Xj = L: 

(20) 



An) 



(n) 



j E {0, . . . , J}, n E Nq. 



^^+l) e („+l) _ s (0)^(n+l) e (n + l) = _ e W V ,W i + ^ s (n+l-*) _ ^ 



(0 



(01 jl (n+1) , / / (n+1) . (n)\ 
- s { Vj7j + 4>J-1 [lj + 7j J 



At the left contact xq = 0, a nonhomogeneous boundary condition can be derived in a 
similar way. Since the potential energy in the left lead is assumed to vanish, the term e^ n ' 
is not needed, and the boundary condition is given by 



(21) 



where 



4 n) + J2 s(n+1 ^ K } -0o7i 



_ o(0) 



(n) 

7o : = ex P 



-2m arctan 



At 
2fr 



-E 



n e Nn. 



We summarize: The Crank-Nicolson scheme (14) with the nonhomogeneous discrete 
transparent boundary conditions (20) and (21) reads as 

(22) B^ n+1) = C^ {n) + S n \ 

where ijj^ = (V>o , ■ ■ ■ > ^j^Y ■, d — (^o > 0, . . . , 0, c^) 1 ". Furthermore, I? is a tridiagonal 
matrix with main diagonal (-s (0) , ii2 - 2 + iuV r 1 (n+1/2) , . . . , - 2 + wVj"^, -s(°>e( n+1 >), 
upper diagonal (1, . . . , 1), and lower diagonal (1, . . . , l,e^ n+1 ^); C is a tridiagonal matrix 



14 



JAN-FREDERIK MENNEMANN, ANSGAR JUNGEL, AND HANS KOSINA 



with main diagonal (0,iR + 2 — wVi n+1 ^ 2 \ . . . ,iR + 2 — wVj 1 ^ 1 ' , 0) , upper diagonal 
(— 1, . . . , —1) and lower diagonal (—1, . . . , —1, — e^ n ')\ furthermore, 

n 

(23) 4 n) = E s(n+1 ~ e) K - <fc>7$°) - s^oli n+1) + <Pi (7? +1) + 7?°) , 

n 

(24) df = £ S <«+M> (^cW - - ^°V,7i n+1) + 0,-1 (7i n+1) + 7?°) • 

3.2. Fast evaluation of the discrete convolution terms. In the subsequent simula- 



tions, scheme (22) has to be solved in each time step and for every wave function ip = ipk, 
k G /C. We recall that the kernel coefficients s( n > need to be calculated only once as they 
do not depend on the wave number k. Let N denote the number of time steps. For each 
k G /C, we require order O(N) storage units and 0(N 2 ) work units to compute the discrete 



convolutions in (23) and (24). For this reason, simulations with several ten thousands of 
time steps are not feasible. To overcome this problem, one may truncate the convolutions 
at some index, since the decay rate of the convolution coefficients is of order 0(n~ 3 / 2 ) [16j 
Section 3.3]. The drawback of this approach is that still more than thousand convolution 
terms are necessary to avoid unphysical reflections at the boundaries. 

The problem has been overcome in [6] by approximating the original convolution coef- 
ficients and calculating the approximated convolutions by recursion. More precisely, 
approximate s^ n ' by 



s( n \ n < v, 



such that 



(25) C< n >(it) := ^^'^ ~ E s(n " 



i=\ i=\ 



can be evaluated by a recurrence formula which reduces the numerical effort drastically. 
As in [6j, we set v = 2 to exclude and from the approximation. In fact, does 
not appear in the original convolutions, whereas is excluded to increase the accuracy. 

Let A G N. The set {b , q . . . , b\, q\} is computed as follows. First, define the formal 
power series 

h ( x ) := S M + s ^x + s^x 2 + ■■■ + s ("+2A~iV+2A-i + . . . , \ x \ < i. 

The first (at least 2A) coefficients are required to calculate the [A— 1|A]-Pade approximation 
of h, h(x) := Pa-i(x)/Q\(x), where Pa-i and Q\ are polynomials of degree A — 1 and 
A, respectively. If this approximation exists, we can compute its Taylor series h(x) = 
lf> v > + ?^ +1 )a; + • • • , and by definition of the Pade approximation, it holds that 

s (n) = s (n) f or a n n G {u,u + 1, ... ,u + 2A — 1}. 
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It can be shown that, if Q\ has A simple roots qi with |g^| > 1 for all i G {1, . . . , A}, the 
approximated coefficients are given by 

(26) 5 (n) = E^r", bf.= - P ^^-^ :0 ' n>u,te{l,...,A}. 

Summarizing, one first computes the exact coefficients s^°\ . . . , s( y + 2A_1 ) followed by the 
[A — 1|A]-Pade approximation. Then one determines the roots of Q\, yielding the numbers 



qi, . . . , <7a. Finally, one evaluates (26) to find the coefficients bo, ... , 6a. We stress the fact 
that these calculations have to be performed with high precision (2A — 1 mantissa length) 
since otherwise the Pade approximation may fail (see [B]). We employ the Python library 
mpmath for arbitrary-precision floating-point arithmetics [23]. As an alternative, one may 
use the Maple script from [BJ Appendix]. 

A particular feature of this approximation is that it can be calculated by recursion. More 



precisely, for n > v + 1, the function C^ n '(u) in (25) can be written as 

CW(«) = X>W(«), 
1=1 

with 

C < ( l \u)=qJ l C l t- 1 \u) + hqfu^- 2 \ n>v + l, C { e u \u) = 0. 



Hence, the discrete convolutions in (23) and (24) are approximated for n > v = 2 by 



(27) J2 s(n+1 ~ e)u(e) * iCin+1) ( u ) + s 



<D U (») 



i=i 



whereas the exact expressions are used for n = and n — 1. As a result, the storage for 
the implementation of the discrete transparent boundary conditions reduces from O(N) to 
0(A). Even more importantly, the work is of order O(AN) instead of 0(N 2 ). 

Obviously, the quality of the approximation depends on A. By construction, we have 
s (n) _ j(n) £ Qr a jj n g {0, . . . , 2A + z/ — 1} but approximates very well even if n is 



much larger [6]. We illustrate in Section 4.3 that the convergence of the complete transient 



algorithm with respect to A is exponential. 

3.3. The complete transient algorithm. In the previous sections, we have explained 
the approximation of the transient Schrodinger equation with discrete transparent bound- 
ary conditions for given potential energy V = Vb ar r + Keif- Here, we make explicit the 
coupling procedure with the Poisson equation for the selfconsistent potential 

d 2 V Rc}f e 2 



Mi 'n[V self ] - n D ), x e (0,L), V sel{ (0,t) = 0, ^ elf (L,t) = -eU{t), 



dx 2 

with the electron density 



n[V se if](x,t) = / g{k)\%p k (x 1 t)\ 2 dk. 
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According to the Crank-Nicolson scheme, a natural approach would be to employ a two- 
step predictor-corrector scheme. More precisely, let {ipk}keK {4>k }keK be propagated 



for one time step using to obtain v}*^. Then one uses V^ elf 



(n+l/2) 



UK 



(n) 



r(*)^ 



propagate {ip^}k&/c {i'T^^^keic again. This procedure doubles the numerical effort and 
is computationally too costly. As an alternative, the scheme V 



(n+l) 



2 V ' self + Keif , 



to 



(n+l/2) 
self 



2V 



(a) 



V, 



(n-1/2) 



can be employed (as in [31]). We found in our simulations that the most simple approach, 
Klvf 1- := Kdf ' gi yes essentially the same results as the above schemes. The reason is 



that the electron density evolves very slowly compared to the small time step size which is 
needed to resolve the fast oscillations of the wave functions. Hence, the variations of V^ e if 
are small. Similarly, the right boundary condition of the Poisson equation can be replaced 
by Kseif(L) = —eU(nAt) if the applied voltage varies slowly. This is used in the circuit 
simulations of Section [Sj 

The complete transient algorithm is presented in Figure [5j 



input: V B 



V sel f {{(f)k}k€ic) 



in }k€K 



{<t>k}h 



tic 



v. 



(n) 



V. 



(n+l) 



M n+1) W <- WfW 



n[Keif]^AfcE t6K ffWl4" +1) l 2 



d 2 T/(»+l) _ e 2 



dx' 1 self 



(n[V aEl{ ) - n D ) , V^ 11 (0) = , (L) = ~eU^ 




Figure 5. Flow chart of the transient scheme. 



3.4. Discretization parameters. We choose K = 3000 for the number of wave functions 
as in the stationary simulations and At = lfs (fs = femtosecond) for the time step size. 
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With the maximal kinetic energy of injected electrons hiou = h 2 k 2 M /(2m), where ku is 
the maximal wave number, the period is computed according to r M = In juj M . Thus, the 
fastest wave oscillation is resolved by TM/At ps 18.5 time steps. The space grid size is 
chosen to be Ax = 0.1 nm. Consequently, the smallest wave length \m = 27r jku ~ 10 nm 
is resolved by approximately 100 spatial grid points. Furthermore, we take A = 70 for 
the approximation parameter of the discrete convolution terms. This choice results from 



a numerical convergence study presented in Section 4.3 



It is important to note that the wave functions which are propagated using the fast 



evaluation of the approximated discrete convolution terms (27) practically coincide with the 



wave functions which are propagated using the exact convolutions (23)-(24) (see Section 



4.3). Employing the exact convolutions, however, is equivalent to solving the Crank- 
Nicolson finite difference equations of the whole space problem. Considering that the 
electron density evolves smoothly in space and time, it is clear that the error of the complete 



transient algorithm (see Section 3.3[ ) is determined by the Crank-Nicolson finite difference 



scheme. A global error estimate, together with a meshing strategy depending on a possibly 
scaled Planck constant h is given in [TJ. The calculations in this article are performed in 
SI units without any scaling. 

3.5. Details of the implementation. The final solver is implemented in the C++ pro- 
gramming language using the matrix library Eigen [22] for concise and efficient compu- 
tations. As we are interested in simulations with a very large number of time steps A" 
(e.g., A" = 100 000), some sort of parallelization is indispensable. We employ the library 
pthreads to realize multiple threads on multi-core processors with shared memory. The 



most time consuming part in the transient algorithm (see Section 3.3) is the propagation 
of the wave functions and the calculation of the electron density. Since the wave functions 
evolve independently of each other, this task can be easily parallelized. At every time step, 
we create a certain number of threads (usually, this number equals the number of cores 
available). To each thread, we assign a subset of wave functions which are propagated as 
described above. Before the threads are joined again, each thread computes its part of the 
electron density. All these parts provide the total eletron density which is used to solve 
the Poisson equation in serial mode. The simulations presented below have been carried 
out on an Intel Core 2 Quad CPU Q9950 with 4 x 2.8 GHz. 



4. Numerical experiments 

We present three numerical examples. The first example shows the importance to provide 
a complete compatible discretization of the open Schrodinger-Poisson system. The second 
numerical test shows the time-dependent behavior of a resonant tunneling diode, which 
allows us to identify three physical regions. In the third experiment, we investigate the 
convergence of our solver with respect to the parameter A which appears in the context of 
the fast evaluation of the discrete convolution terms. 

4.1. First experiment: Constant applied voltage. We compute the stationary solu- 
tion to the Schrodinger-Poisson system with an applied voltage of U = 250 mV. At this 
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voltage, the current density achieves its first local maximum. We apply the transient algo- 
rithm of Section [3] until t = 25 fs, keeping the applied voltage constant. Accordingly, the 
stationary solution should be preserved and the current density J con d, defined in (12), is 
expected to be spatially constant. 

The ad-hoc discretization ( |T9| ) is employed using the time step sizes At = lfs, 0.5 fs, 
0.25 fs. We observe in Figure [6] that the current density is not constant. The reason is 



that the discretization (19) is not compatible with the underlying finite-difference scheme. 



The distortions are reduced for very small time step sizes but this leads to computationally 



expensive algorithms. In contrast, with the discrete gauge change of Section 3.1 the current 
density is perfectly constant even for the rather large time step At = 1 fs; see Figure [6) 
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CD 
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dt = 1 fs 
dt = 0.5 fs 
■ dt = 0.25 fs 
• dt = 1 fs 
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Figure 6. Conduction current density in a resonant tunneling diode at 
t = 25 fs for a constant applied voltage of U = 250 mV. Discretizations 



using the ad-hoc discretization (19) of the analytical boundary conditions 
yield strongly distorted numerical solutions (broken lines). In contrast, the 
conduction current density computed with our solver is perfectly constant 
(solid line). 



We mention that the transient solution is also distorted if the scattering states as initial 
wave functions are computed from an ad-hoc discretization of the continuous boundary 
conditions (|5j and For stationary computations, spurious reflections due to an in- 
consistent discretization play a minor role but they become a major issue for transient 
simulations. 

4.2. Second experiment: Time-dependent applied voltage. For the second numer- 
ical experiment, we consider a time-dependent applied voltage. The conduction current 
density is no longer constant but the total current density is expected to be conserved. 
We recall that the total current density J to t = </ CO nd + dD /dt is the sum of the conduction 
current density J con d and the displacement current density dD/dt. Here D denotes the 
electric displacement field which is related to the electric field E by D = tQe r E. Indeed, 
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replacing the electric field by the negative gradient of the potential we obtain 

3D e e r d 



dt " e dt 



self- 



The temporal and spatial derivatives are approximated using centered finite differences. 
Ampere's circuital law V x H = J tot for the magnetic field strength H yields 

div J tot = div(V x H) =0, 

and hence, in one space dimension, J tot is constant in space. 

The following simulation demonstrates that the total current density is a conserved 
quantity in the discrete system as well. First, we compute the equilibrium state using an 
applied voltage of U = V. This solution is then propagated using a raised cosine function 
for the applied voltage 

t/(t) = y (l - cos ^ , 0<t<l P s, 

where Uq = 0.25 V and T = 2ps. At later times, t > lps, U(t) = Uq is kept constant. 
Conduction, displacement, and total current density at different times are depicted in the 
left column of Figure [7} As can be seen, the total current density is perfectly conserved 
at all considered times. The change of the charge density dp/dt is illustrated in the right 
column of Figure [7| In our model, p is given by p = e (n D — n). 

The time-dependence of the total current density in response to the applied voltage is 
shown in Fig. [8j We can identify three different regions in the temporal behaviour, each 
of which is governed by a different physical mechanism. 

Region I: Capacitive behavior. When the applied voltage increases during the first picosec- 
ond, the resonant tunneling diode behaves mainly like a parallel plate capacitor. This can 
be clearly seen in the top left panel of Figure [7} In the region of the double barrier, the 
displacement current gives the dominant contribution to the total current, whereas the 
conduction current is small. The top right panel of Figure [7] shows a build-up of negative 
charge before the left barrier and of positive charge after the right barrier. The formation 
of opposite charges on the two sides of the double barrier results in the formation of an 
electric field between the two regions of opposite charge density. This field is necessary to 
accomodate the externally applied voltage. Figure [8] shows that the current closely follows 
the time derivative of the applied voltage: 

dU 7iCU . (2%t\ 

Jcond ~ C — — = — — Sin — — . 

dt T \T J 

This expression allows us to estimate the apparent capacitance C. The maximum current 
density occurring at t = T/4 = 0.5 ps takes approximately the value 1.2 • 10 9 Am -2 . We 
compute C = TJ/tiUq = 3.06 • 10 -3 Fm -2 . Equating this value to the parallel plate 
capacitance, C = e e r /d, we find the average separation of the opposite charge densities 
to be d — 33.1 nm. 
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Figure 7. Left column: Total current density J tot = J con d + dD/dt, con- 
duction current density J := J COQ d, and displacement current density dD/dt 
versus position at different times. Right column: Temporal variation dp/dt 
of the charge density versus position. 



Region II: Plasma oscillations. During the second picosecond, a strongly damped oscilla- 
tion occurs in the current density. From Figure [8j we estimate five oscillations to occur 
during one picosecond, which relates to a period of about 200 fs. It is believed that these 
are plasma oscillations which were excited by the rapidly changing applied voltage U. As 
soon as the transient phase of U(t) is over and U(t) is kept constant at Uq for t > lps, the 
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Figure 8. Applied voltage and total current density versus time in different 
scalings. 



excitation vanishes and the oscillations fade out quickly. As a rough estimate we calculate 
the plasma frequency u p for a classical electron system of uniform density: 

2 

U = • 

y ms 

Note that in the resonant tunneling diode the density is neither uniform nor is it governed 
by the classical equations of motion. Nevertheless, we may use this expression to estimate 
the order of magnitude of the time constant associated with this effect. Since plasma 
oscillations usually occur in the high-density regions of a device, we set n = n l D = 10 24 m -3 
and obtain r p = 2ir/u p = 111.4 fs. This value is of the same order as the 200 ps estimated 
above, which is a strong indication that the physical effect observed here is a plasma 
oscillation. 

Region III: Charging of the quantum well. For t > 2 ps, an exponential increase in the 
current can be clearly observed in Figure |8j Below 2 ps we see a superposition of both the 
exponential current increase and the plasma oscillations. The origin of this effect can be 
understood from the right panels of Figure [7j Negative charge builds up in the quantum 
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Figure 9. Number of electrons in the quantum well versus time. In Region 
III (t > 2ps) this number clearly follows an exponential law. U a2 , a5 denotes 
the temporal variation of the voltage between x = a<i and x = a§. 



well. This charge results from electrons tunneling through the left barrier into the quantum 
well. In this context, we note that the temporal variation of the voltage between the left 
and right end points a 2 and as of the double-barrier structure, respectively, follows closely 
the variation of the applied voltage U and hence, it is practically constant for t > 1 ps 
(see Figure [9j). The rate \dp/dt\ decreases with time as can be seen by the snapshots at 
t — 1.5 ps and t = 3 ps. We calculate the number of electrons residing in the quantum well: 

ra-b 

N(t) := / n(x,t)dx. 

Since the charging process is expected to show an exponential time dependence, we assume 
the following exponential law for N(t) and extract the free parameters r and iV^: 

N(t) =N oa + (Nih) - AU e- {t -^l\ 

In Figure [9j the difference \N(t) — N^] is plotted, which decays to zero with an extracted 
time constant of r = 1.25 ps. 

This time scale is related to the life time of a quasi-bound state. At U = 0.25 V, the 
current- volt age characteristic has its first maximum, which means that the first resonant 
state in the quantum well is carrying the current. The life time of this resonant state can 
be extracted from the width of the resonance peak in the transmission coefficient. Figure 



10 depicts the transmission coefficient of the double-barrier structure at U — 0.25 V. The 
transmission coefficient is defined as the ratio between the transmitted and the incident 
probability current density jtrans and ji nc . In terms of the amplitude and the wavenumber 
of the transmitted and the incident wave, it reads: 



|jtrans 




A 

^Hrans 


2 h 
"'trans 


|jinc| 
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Figure 10. Transmission coefficient of the double-barrier structure at U 
0.25 V in different scalings. 



Extracting AE, the half width at half maximum of the first transmission peak, the life 
time of the resonant state can be estimated as follows [25] : 

h 



T 



2AE 

At U = 0.25 V we find 2AE = 5.31 • 10~ 4 eV and thus r = 1.24 ps. This value is very close 
to the time constant of r = 1.25 ps extracted from the exponential charge increase in the 
quantum well, which is the cause for the observed exponential current increase. 

4.3. Third experiment: Convergence in A. Finally, we study the convergence of the 



complete transient algorithm detailed in Section 3J3 with respect to the parameter A which 
appears in the context of the fast evaluation of the discrete convolution terms. For this 
purpose, we repeat the last experiment with different values of A. We compare the results 
with those obtained from the algorithm which uses the discrete transparent boundary 



conditions with the exact convolutions (23)-(24). Since the computation of the reference 
solution is extremely expensive, we restrict the experiment to the final time t = 1.5 ps. 
The conduction current densities at t — 1.5 ps for two different values of A and for the 

? 2 -norm for 



11 



(left). The relative error in the 
right). We observe that the relative error 



reference solution are depicted in Figure 
increasing values of A is shown in Figure 11 
decreases exponentially fast. Thus, using a relatively small value of A practically yields 
the same results (at dramatically reduced numerical costs) as if the discrete transparent 
boundary conditions with the exact convolutions were used. 

5. Circuit simulations 

In this section, we simulate a high-frequency oscillator consisting of a voltage source U e , 
a resistor with resistance R, an inductor with inductance L, a capacitor with capacity C, 



and a resonant tunneling diode RTD; see Figure [12j Each element of the circuit yields one 
current-voltage relationship, 

(28) U R = RI R , U L = Li L , I C = CU C , J RTD = /(C/ R td). 
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Figure 11. Conduction current density at t 
terror for increasing A. 



1.5 ps (left) and relative 



The last expression is to be understood as follows. Given the applied voltage £/rtd at 
the tunneling diode, the current /rtd(^) = AJ tot (t) is computed from the solution of the 
time- dependent Schrodinger-Poisson system. Here, A = 10 -11 m 2 is the cross sectional 
area of the diode and J tot is the total current density. In the simulations we use R = 5 Q, 
L = 50pH, and C = lOfF. 

7? 




Figure 12. High-frequency oscillator containing the resonant tunneling 
diode RTD. 



According to the Kirchhoff circuit laws, we have 



(29) 



U e = U R + U R TD + U L , U } 



RTD 



u, 



c, 



Ir, 



Combining (28) and (29), we find that 

CUbtd = CUc = Ic = II — ^rtd, 
LI L = U L = U e -U R - C/rtd = 



'RTD 



+ Ic- 



U e - RI R - U RTD = U e - RI L - U } 



RTD- 
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Consequently, we obtain a system of two coupled ordinary differential equations, 
(30) ±(U^\ = (0 h^flk^ + (-il^ 



it \ il J ~ \-i -fj v J i y v m) j ■ 

The time-step size At is very small compared to the time scale of the variation of the 
potential energy and the variation of the current flowing through the diode. Hence, using 



the same time step for the time integration of (30), we can resort to an explicit time- 
stepping method. We choose the simplest one, the explicit Euler method. Alternatively, 
one may employ an implicit method, but we observed that both methods yield essentially 
the same results. 

First circuit simulation. In the first simulation, the RTD solver is initialized with the steady 
state corresponding to £/rtd(£) = for all t < 0. The external voltage U e is assumed to 



be zero for t < 0, and the initial conditions for (30) are £/rtd(0) = and iz,(0) = 



For t G [10, 20] ps, the external voltage is increased smoothly to 0.275 V and then kept 



constant (see Figure 13). This value is between the voltages where the stationary current 
density reaches its local maximum and minimum (see Figure |4|). The time evolution of the 



voltage and the current at the RTD are depicted in Figure 13 It is clearly visible that 



the system starts to oscillate. Furthermore, the potential energy, electron density, current 



densities, and the temporal variation of the total charge dp/dt are shown in Figure 14 for 
four different times from the interval \pi = 77.7,t 2 = 87.2] ps, which covers exactly one 
oscillation. Around 2 ps after the beginning of the period, the electron density within the 
quantum well in [65, 70] nm becomes minimal (first row). After some time, we observe a 
build-up of negative charge in the quantum well with dp/dt < (second row). At about 
t = 84.6 ps the electron density reaches its maximum value (third row). Subsequently, the 
electrons leave the quantum well again and dp/dt > in [65, 70] nm (fourth row). The 
frequency of the oscillations is approximately 105 GHz which corresponds qualitatively 
to frequencies observed in standard double-barrier tunneling diodes |T4] . The temporal 
evolution of the physical quantities in the circuit is animated in the video available at 
http : //www. asc . tuwien. ac . at/~mennemann/projects .html. 

Second circuit simulation. In this experiment, the external voltage U e is kept fixed for all 
times. At times t < 0, the circuit contains the voltage source, resistor, and RTD only. We 
initialize the transient Schrodinger-Poisson solver with the steady state corresponding to 
^rtd(^) = 0.275 V for all t < 0. To compensate for the voltage drop at the resistor, the 
external voltage is set to 

U e (t) = Rljxroit) + U KTD (t), t < 0. 

At time t = 0, the capacitor and inductor are added to the circuit. In order to avoid 
discontinuities in the voltages, we charge the capacitor with the same voltage wich is 
applied at the RTD before the switching takes place, Uc(t) = £7rtd(£) for t < 0. For 
similar reasons, we set the current flowing through the inductor to the current flowing 
through the RTD, I Lit) = Jrtd(^) for t < 0. This configuration corresponds to the 
equilibrium state. Therefore, one would expect that the system remains in its initial state 
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25 50 75 100 

time in ps 

Figure 13. First circuit simulation: Voltage C/rtd an d current Irtd 
through the resonant tunneling diode versus time. 



for all time. However, the equilibrium is unstable and a small perturbation will drive the 
system out of equilibrium. In fact, numerical inaccuracies suffice to start the oscillator. 
However, we accelerate the transient phase by perturbing by the value 5 • 1CT 6 A for 



t < 0. The numerical result is presented in Figure 15 The simulation took less than 4 



hours computing time on an Intel Core 2 Quad Core Q9950 with 4 x 2.8 GHz. 
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Figure 14. First circuit simulation: Electron density, potential energy, cur- 
rent densities, and variation of the electron density versus position in the 
RTD at four different times. 
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Figure 15. Second circuit simulation: Applied voltage and current at the 
resonant tunneling diode versus time. 
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